Logistic Map 综合

2025-10-8

课程代码均在如下 GitHub 仓库开源：https://github.com/GHe0000/scientific-programming-intro

注：由于计算和图片较多，从头运行大概需要 40 秒.

目录

• 作业要求
• 解析求解周期一和周期二轨道及其稳定的参数范围
  • 基本思路
  • 周期一
  • 周期二
  • 为何会出现 Overflow
• 全周期分叉图
  • 分叉图
  • 周期三
  • 周期五
• 全参数范围的李雅普诺夫指数
• 倍周期分叉点、$r_{\infty}$、费根鲍姆常数
  • 倍周期分叉点
  • $r_{\infty}$
  • 验证倍周期分叉点和 $r_{\infty}$
  • 费根鲍姆常数
• 选做内容：数值计算李雅普诺夫指数
• 结论
• 参考资料
• 致谢
• 附一：问题
• 附二：学习心得

作业要求

对映射 $x_{n+1} = r - x_n^2,\quad r\in[0,2],\quad x\in[-2,2]$ 做如下分析：

• 解析求解周期一和周期二轨道及其稳定的参数范围，并用数值验证
• 画出全参数分叉图，并放大显示某个周期三、周期五区域
• 利用解析公式计算全参数范围的李雅普诺夫指数，是否跟分岔图解释的现象一致？
• 利用全参数范围的李雅普诺夫指数确定倍周期分岔点以及 $r_{\infty}$，利用倍周期分岔点的参数计算费根鲍姆常数 $\delta$，是否跟 Logistic map 一致？

选做内容（只选其一）：

• 定量分析临界慢化: 倍周期分岔过程中，当参数趋向分岔点时，轨道趋向吸引子的速度会变慢，定量数值研究收敛速度与参数到分岔点距离的依赖关系，不同的分岔点规律是否一致？
• 数值计算李雅普诺夫指数：根据课堂介绍的数值方法计算李指数，并挑选若干参数跟解析结果对比。

要求：排版简洁规范，结构完整。要有必要的文字描述解题思路、编程逻辑、现象分析和总结。

解析求解周期一和周期二轨道及其稳定的参数范围

基本思路

对于映射 $x_{n+1} = r-x^2_{n}$，我们定义迭代函数：

$$ f(x) = r - x^2,\quad x_{n+1} = f(x_n) $$

同时我们定义：

$$ f^n (x) = \underbrace{\left( f \circ f \circ \dots \circ f \right)}_{n \text{个} f} (x)$$

要求解周期一的轨道，就是要求解 $f(x)$ 的不动点，即求解方程 $x^* = f(x^*)$，而不动点的稳定性可以通过局部稳定判据进行判断，也就是计算 $|f'(x^*)|$ 并和 $1$ 进行对比，求解出 $|f'(x^*)|<1$ 的 $r$ 的范围就是周期一轨道的稳定的参数范围.

而周期二的轨道，就是要求解 $f^2(x)$ 的不动点，即求解方程 $x^* = f^2(x^*)$，并类似地计算 $|(f^2)' (x^*)|$ 进行局部稳定性的判断.

周期一

解析求解

这里我们先导入需要的库，并针对 Jupyter Notebook 进行一些设置.

import numpy as np # 数值计算库
import numba as nb # 引入 jit 来加速函数（如果需要）
import sympy as sym # 符号计算库
import matplotlib as mpl
import matplotlib.pyplot as plt # 图像绘制
from matplotlib.backends.backend_agg import FigureCanvasAgg # 设置图像后端用
import ipywidgets as ipw # 交互控件
from IPython.display import display, Math, Latex # 打印数学公式

# 设置随机数种子确保结果的可复现性
np.random.seed(3407)

# 使用 mathjax 来在 jupyter notebook 显示数学公式
sym.init_printing(use_latex='mathjax')

# 设置 matplotlib 绘制的图像嵌入到 jupyter notebook 的方式
%matplotlib widget

# 一个工具函数，可以让静态图片不经过 widget 直接嵌入 Jupyter notebook
# 这样图片可以直接存到 nb 文件里，而不是需要运行才能显示（类似 inline）
def display_inline(fig):
    fig.set_canvas(FigureCanvasAgg(fig))
    display(fig)
    plt.close(fig) # 释放 fig，减小资源消耗

# 一个工具函数，用于在 Jupyter 中快速实现公式和字符串混合输出
def display_math(*args):
    parts = []
    for a in args:
        if isinstance(a, str):
            parts.append(a)
        else:
            parts.append(sym.latex(a))
    display(Math("".join(parts)))

# 设置 matplotlib 使用的字体，避免出现中文问题
plt.rcParams['font.family'] = ['SimSun', 'Times New Roman']
plt.rcParams['mathtext.fontset'] = 'cm'
# 设置默认坐标轴字体大小
plt.rcParams['axes.labelsize'] = 14

这里我们定义一个函数 fn，从而方便地得到任意的 $f^n (x)$ 函数. 并以此定义 $f(x)$（记为 $f^1(x)$）

def fn(r, x, n):
    f = lambda r,x: r-x**2
    for _ in range(n):
        x= f(r,x)
    return x

sym.var("r,x") # 定义符号
f1 = fn(r,x,1)
display_math(r"f(x)=",f1)

$\displaystyle f(x)=r - x^{2}$

计算 $f(x)$ 的不动点 $x^*$，并计算 $f'(x)$，将 $x^*$ 带入 $f'(x)$：

sol = sym.solve(f1-x,x)
df1 = sym.diff(f1,x)

display_math(r"\text{函数 } f(x)=",f1,f"\\text{{ 一共找到 {len(sol)} 个不动点：}}")
args_list = []
for s in sol:
    args_list.extend([r"x^*=", s, r",\quad "])
display_math(*args_list)

display_math(r"f(x) \text{ 的导数为 }f'(x)=",df1,r"\text{, 带入 } x^* \text{ 有：}")
for s in sol:
    display_math(r"f'\left(", s, r"\right)=", df1.subs(x,s))

$\displaystyle \text{函数 } f(x)=r - x^{2}\text{ 一共找到 2 个不动点：}$

$\displaystyle x^*=- \frac{\sqrt{4 r + 1}}{2} - \frac{1}{2},\quad x^*=\frac{\sqrt{4 r + 1}}{2} - \frac{1}{2},\quad $

$\displaystyle f(x) \text{ 的导数为 }f'(x)=- 2 x\text{, 带入 } x^* \text{ 有：}$

$\displaystyle f'\left(- \frac{\sqrt{4 r + 1}}{2} - \frac{1}{2}\right)=\sqrt{4 r + 1} + 1$

$\displaystyle f'\left(\frac{\sqrt{4 r + 1}}{2} - \frac{1}{2}\right)=1 - \sqrt{4 r + 1}$

由于 $r\in[0,2]$，因此 $\sqrt{4r+1}+1>1$，不动点：

$$ x^* = \frac{-\sqrt{4r+1}-1}{2} $$

一定是不稳定的. 而不动点：

$$ x^* = \frac{\sqrt{4r+1}-1}{2} $$

若是稳定的，则需要：

$$ \left| 1-\sqrt{4r+1} \right| < 1 $$

解得：$ r < \frac{3}{4} $. 因此，当 $0<r<\frac{3}{4}$ 时，系统有稳定的周期一轨道：

$$ \dots \to \frac{\sqrt{4r+1}-1}{2} \to \frac{\sqrt{4r+1}-1}{2} \to \dots $$

数值验证

为了方便后面绘制迭代过程，我们将绘制迭代过程抽象成一个函数：

def plot_orbit(ax, map_func, x0, r, n, **kwarg):
    n_arr = np.arange(n+1)
    x_arr = np.zeros((n+1, len(x0)), dtype=np.float64)
    x_arr[0] = x0.copy()
    for i in range(n):
        x_arr[i+1] = map_func(r, x_arr[i])
    for i, x in enumerate(x0):
        ax.plot(n_arr, x_arr[:,i], label=f"$x_0$={x:.3f}",  **kwarg)
    return ax

给定 $r=0.5$，我们在 $x\in[-2,2]$ 之间均匀采样 $5$ 个点，绘制其迭代过程：

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 0.5
ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(-2,2,5), r, n, markersize=4, marker='o', linestyle='--')
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n})")
ax.axhline((np.sqrt(4*r+1)-1)/2, color='r', linestyle=':', label=r"预期值：$(\sqrt{4r+1}-1)/2$")
ax.grid()
ax.legend()

display_inline(fig)

/tmp/ipykernel_585326/3616708647.py:5: RuntimeWarning: overflow encountered in square
  ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(-2,2,5), r, n, markersize=4, marker='o', linestyle='--')

可以发现，出现了 Overflow 的警告，同时图像的坐标轴被放很大，我们设置一下人为忽略 Overflow 的警告（在后面的章节会对其进行具体分析），并限制一下 $y$ 轴的范围：

为了显示出现了 Overflow，上面代码保留了出现的警告，后文中会经过处理来消除.

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 0.5
with np.errstate(over='ignore', invalid='ignore'):
    ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(-2,2,5), r, n, markersize=4, marker='o', linestyle='--')
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n})")
ax.axhline((np.sqrt(4*r+1)-1)/2, color='r', linestyle=':', label=r"预期值：$(\sqrt{4r+1}-1)/2$")
ax.grid()
ax.legend()

display_inline(fig)

我们选取更多的 $x_0$ 查看变化过程：

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 0.5
with np.errstate(over='ignore', invalid='ignore'):
    ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(-2,2,20), r, n, markersize=4, marker='o', linestyle='--')
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，20 个初始值)")
ax.axhline((np.sqrt(4*r+1)-1)/2, color='r', linestyle=':', label=r"预期值：$(\sqrt{4r+1}-1)/2$")
ax.grid()
#ax.legend()

display_inline(fig)

可以看到，对于某些 $x_0$，$x_n$ 最终到达稳定的周期一轨道，而对于另外的 $x$，$x_n$ 最终会到达 $- \infty$

我们根据 $x_n$ 最后是否收敛到周期一轨道对线条进行染色，并加大 $x_0$ 采样的密度：

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 0.5
x0 = np.linspace(-2,2,50)
map_func = lambda r,x: r-x**2
with np.errstate(over='ignore', invalid='ignore'):
    n_arr = np.arange(n+1)
    x_arr = np.zeros((n+1, len(x0)), dtype=np.float64)
    x_arr[0] = x0.copy()
    for i in range(n):
        x_arr[i+1] = map_func(r, x_arr[i])
    for i, x in enumerate(x0):
        line_color = "red" if np.isinf(x_arr[-1,i]) else "green"
        ax.plot(n_arr, x_arr[:,i], label=f"$x_0$={x:.3f}",  markersize=4, marker='o', linestyle='--', color = line_color)
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，50 个初始值，绿线被吸引)")
ax.axhline((np.sqrt(4*r+1)-1)/2, color='k', linestyle=':', label=r"预期值：$(\sqrt{4r+1}-1)/2$")
ax.text(9, 0.1, f"稳定的周期一轨道，$\\frac{{\\sqrt{{4r+1}}-1}}{{2}} \\approx {(np.sqrt(4*r+1)-1)/2 :.3f} $", color='k', size=12)
ax.axhline(1.36, color='blue', linestyle=':')
ax.axhline(-1.36, color='blue', linestyle=':')
ax.text(12, 1.43, r"分界线 1, $x_0\approx=1.36$", color='b', size=12)
ax.text(12, -1.33, r"分界线 2, $x_0\approx-1.36$", color='b', size=12)
ax.grid()
#ax.legend()

display_inline(fig)

可以看到，当初始值在两条分界线之间时，$x_n$ 收敛到周期一轨道，反之到达 $-\infty$.

我们将收敛到周期一轨道的 $x_0$ 构成的集合称为该轨道的吸引域，下面我们尝试求解吸引域的范围.

解析估计吸引域

前面我们通过导数对周期一轨道进行稳定性判断的方法依赖于 Taylor 展开，也就是其只能在吸引子附近的邻域内才能使用. 对于 $|f'(x^*)|<1$ 只能说明其在不动点周围的小邻域范围内是吸引的，对于全取值范围内不一定能使用. 因此，要计算（或者估计）吸引域的范围，需要一些新的方法.

首先想到的是数学分析中的 Banach 压缩映射定理，其要求在某个不变的完备度量空间 $(X,d)$，映射 $T$ 是压缩映射，即 $\forall x,y\in X$ 满足：

$$ d(T(x),T(y)) \leq \lambda d(x,y),\quad (0\leq \lambda <1 ) $$

则 $T$ 在 $X$ 上存在且仅存在一个不动点 $x^*$，即 $T(x^*)=x^*$，且 $\forall x_0 \in X$，迭代 $x_{n+1} = T(x_n)$ 均收敛到唯一的不动点 $x^*$，且收敛速度至少是几何级数，即：

$$ d(x_n,x^*) \leq \frac{\lambda^n}{1-\lambda} d(x_1,x_0) $$

---

对于一维系统，取常规的距离度量，给定区间 $I$，由 Lagrange 中值定理，有：

$$ |f(x)-f(y)| = |f'(\xi)| |x-y| \leq \max_{x\in I} |f'(x)| |x-y|,\quad (\xi \in I) $$

因此 Banach 压缩映射条件可以变为：

$$ |f'(x)|\leq \lambda < 1,\ \forall x \in I $$

因此，对于映射 $f(x)=r-x^2$，如果我们能找到一个区间 $I$，使得 $f(I)\subset I$，并且：

$$ |f'(x)| \leq \lambda < 1, \forall x\in I $$

那么根据 Banach 压缩映射定理，所有初始值在 $I$ 的点最终都会收敛到 $I$ 中唯一的不动点. 则区间 $I$ 一定是吸引域的子集.

且显然可以得到推论，如果 $|f'(x)|>1$ 则该点必然是排斥的.

上述过程的麻烦之处在于找到区间 $I$，使得 $f(I)\in I$. 根据对称性，我们假设区间 $I \in (-M,M)$，因此，$M$ 的最大值满足：

$$ r-M_1^2 = M_1,\ r-M_2^2=-M_2,\ M=\max\{|M_1|,|M_2|\} $$

即：

$$ M=\max\left\{ \left| \frac{\pm\sqrt{4r+1}-1}{2} \right| , \left| \frac{\pm\sqrt{4r+1}+1}{2} \right| \right\} = \frac{\sqrt{4r+1}+1}{2}$$

因此，若 $x_0 \in [-M,M]$，能保证在之后的迭代中一定不超过 $[-M,M]$（至于是否收敛到周期一需要计算导数.）

对于 $f(x)=r-x^2$，$|f'(x)|<1$ 解得 $-1/2 < x < 1/2$，同时 $1/2 < M$，因此由 Banach 压缩映射，我们知道 $[-1/2,1/2]$ 是周期一轨道的吸引域的子集.

这个估计非常保守，这是因为此范围保证了 $f(x)$ 是压缩映射，但是实际上有可能 $f(x)$ 不是压缩映射，但 $f^n(x)$ 是压缩映射，显然这种情况最后也会收敛到不动点，因此使得估计的范围过小.

而 $[-M,M]$ 保证了一旦 $x_n$ 落入这个范围，则 $x_n$ 之后一定保持在这个范围内，根据前面的数值验证，对于存在稳定的周期一轨道的情况，如果 $x_n$ 不趋向 $\pm\infty$，则 $x_n$ 最终收敛于周期一轨道，因此我们可以将 $[-M,M]$ 作为周期一轨道的一个估计.

（但证明 $\forall x_0 \in [-M,M]$ 均收敛到周期一轨道需要借助 Banach 定理通过更精细的分析和估计来证明，能力有限不进一步证明.）

我们将估计的吸引域范围绘制在图上.

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 0.5
x0 = np.linspace(-2,2,50)
map_func = lambda r,x: r-x**2
guess = lambda r: (np.sqrt(4*r+1)+1) / 2
with np.errstate(over='ignore', invalid='ignore'):
    n_arr = np.arange(n+1)
    x_arr = np.zeros((n+1, len(x0)), dtype=np.float64)
    x_arr[0] = x0.copy()
    for i in range(n):
        x_arr[i+1] = map_func(r, x_arr[i])
    for i, x in enumerate(x0):
        line_color = "red" if np.isinf(x_arr[-1,i]) else "green"
        ax.plot(n_arr, x_arr[:,i], label=f"$x_0$={x:.3f}",  markersize=4, marker='o', linestyle='--', color = line_color)
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，50 个初始值，绿线被吸引)")
ax.axhline((np.sqrt(4*r+1)-1)/2, color='k', linestyle=':', label=r"预期值：$(\sqrt{4r+1}-1)/2$")
ax.text(9, 0.1, f"稳定的周期一轨道，$\\frac{{\\sqrt{{4r+1}}-1}}{{2}} \\approx {(np.sqrt(4*r+1)-1)/2 :.3f} $", color='k', size=12)
ax.axhline(1.36, color='blue', linestyle=':')
ax.axhline(-1.36, color='blue', linestyle=':')
ax.text(12, 1.43, r"分界线 1, $x_0\approx=1.36$", color='b', size=12)
ax.text(12, -1.33, r"分界线 2, $x_0\approx-1.36$", color='b', size=12)
ax.axhline(guess(r), color='k', linestyle='-.')
ax.text(9, guess(r)-0.25, f"解析估计吸引域上界：$\\frac{{\\sqrt{{4r+1}}+1}}{{2}}\\approx{guess(r):.3f}$", color='k', size=12)
ax.axhline(-guess(r), color='k', linestyle='-.')
ax.text(9, -guess(r)+0.25, f"解析估计吸引域下界：$-\\frac{{\\sqrt{{4r+1}}+1}}{{2}}\\approx{-guess(r):.3f}$", color='k', size=12)
ax.grid()
#ax.legend()

display_inline(fig)

可以发现我们的估计还是比较准确的.

数值计算吸引域

前面我们尝试通过解析的方式估计吸引域，但需要通过精细的分析来估计，非常麻烦且个人能力有限. 目前仅猜测 $[-M,M]$ 是周期一的吸引域. 其中：

$$ M = \frac{\sqrt{4r+1}+1}{2} $$

这里我们将尝试通过数值的方法找到吸引域的范围.

首先，由于 $f(x)$ 是关于 $x=0$ 对称的，根据对称性我们可以合理的认为吸引域也是关于 $x=0$ 对称的. 同时，观察前面的收敛图，可以知道吸引域满足 $(-N,N)$，因此，要计算周期一的吸引域，就是要确定 $N$. 这里使用二分法来寻找.

1. 给定一个初始范围 $(N_1, N_2)$，保证边界 $N$ 落在这个范围内
2. $N_x = (N_1 + N_2) / 2$，以 $N_x$ 为初始值迭代，若 $N_x$ 在周期一轨道或者得到 inf 结束迭代
3. 若在周期一轨道，更新 $N_2 \to N_x$，若得到 inf，更新 $N_1 \to N_x$
4. 回到第二步，直到 $N_2-N_1$ 小于给定范围，则 $N_x$ 为最后结果.

def find_point(check_func, r_arr, tol=1e-6, max_step=50): # 二分查找函数，向量化
    # check_func 1 改变下界, 0 改变上界
    a = r_arr[:, 0].copy()
    b = r_arr[:, 1].copy()
    for _ in range(max_step):
        mid = 0.5 * (a + b)
        check = check_func(mid)
        check = np.asarray(check, dtype=bool)
        a = np.where(check, mid, a)
        b = np.where(check, b, mid)
        if np.all(np.abs(b - a) < tol):
            break
    return 0.5 * (a + b)

def make_check_func(r_arr, tol=1e-3, max_step=100): # 生成判断函数，向量化
    r_arr = np.array(r_arr,ndmin=1)
    x_star = (-1 + np.sqrt(1 + 4*r_arr)) / 2.0
    def check_func(x0):
        x = x0.copy()
        with np.errstate(over='ignore', invalid='ignore'):
            for _ in range(max_step):
                x = r_arr - x**2
        return np.abs(x - x_star) < tol
    return check_func

带入 $r=0.5$，初始边界取 $(0,2)$，计算吸引域边界：

r_arr = np.array([0.5])
r0_guess = np.array([[0,2]])
check_func = make_check_func(r_arr)
n_arr = find_point(check_func, r0_guess, tol=1e-6)
print(f"r=0.5 时正半轴吸引域边界：{n_arr[0]}")

r=0.5 时正半轴吸引域边界：1.366025447845459

可以发现和前面我们通过图像确定的值接近，说明上述算法是有效的，我们在 $r\in[0,3/4]$ 内寻找吸引域边界，并绘图：

n_r = 1000
guess_func = lambda r: (np.sqrt(4*r+1)+1)/2
r_arr = np.linspace(0,3/4,n_r)
r0_guess = np.array([[0,2] * n_r])
check_func = make_check_func(r_arr)
n_arr = find_point(check_func, r0_guess, tol=1e-6)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, n_arr, label="二分数值求解")
ax.plot(r_arr,guess_func(r_arr), "--", label=r"猜测的解析解 $(\sqrt{4r+1}+1)/2$")
ax.grid()
ax.set_xlabel(r"$r$")
ax.set_ylabel(r"吸引域上界 $N$")
ax.set_title("稳定周期一轨道内吸引域上界随 $r$ 变化")
ax.legend()
display_inline(fig)

可以发现前面符合很好，但是到 $0.7$ 附近突然发生了变化，我们挑选 $r=0.74$ 查看情况：

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 0.74
x0 = np.linspace(-2,2,50)
map_func = lambda r,x: r-x**2
guess = lambda r: (np.sqrt(4*r+1)+1) / 2
with np.errstate(over='ignore', invalid='ignore'):
    n_arr = np.arange(n+1)
    x_arr = np.zeros((n+1, len(x0)), dtype=np.float64)
    x_arr[0] = x0.copy()
    for i in range(n):
        x_arr[i+1] = map_func(r, x_arr[i])
    for i, x in enumerate(x0):
        line_color = "red" if np.isinf(x_arr[-1,i]) else "green"
        ax.plot(n_arr, x_arr[:,i], label=f"$x_0$={x:.3f}",  markersize=4, marker='o', linestyle='--', color = line_color)
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，50 个初始值，绿线被吸引)")
ax.axhline((np.sqrt(4*r+1)-1)/2, color='k', linestyle=':', label=r"预期值：$(\sqrt{4r+1}-1)/2$")
ax.text(9, 0.1, f"稳定的周期一轨道，$\\frac{{\\sqrt{{4r+1}}-1}}{{2}} \\approx {(np.sqrt(4*r+1)-1)/2 :.3f} $", color='k', size=12)
ax.axhline(guess(r), color='k', linestyle='-.')
ax.text(9, guess(r)-0.25, f"解析估计吸引域上界：$\\frac{{\\sqrt{{4r+1}}+1}}{{2}}\\approx{guess(r):.3f}$", color='k', size=12)
ax.axhline(-guess(r), color='k', linestyle='-.')
ax.text(9, -guess(r)+0.25, f"解析估计吸引域下界：$-\\frac{{\\sqrt{{4r+1}}+1}}{{2}}\\approx{-guess(r):.3f}$", color='k', size=12)
ax.grid()
#ax.legend()

display_inline(fig)

可以发现我们的估计和实际的分界还是吻合的. 同时当 $r$ 接近分叉点 $3/4$ 出现了临界慢化现象，导致二分查找中判断是否收敛时实际是收敛的，但是因为临界慢化现象，结束计算时仍未收敛到判定为收敛的小范围内，从而导致二分查找出现问题.

这里我们加大收敛的判定范围，同时增大判断时的迭代过程重新绘图.

n_r = 1000
guess_func = lambda r: (np.sqrt(4*r+1)+1)/2
r_arr = np.linspace(0,3/4,n_r)
r0_guess = np.array([[0,2] * n_r])
check_func = make_check_func(r_arr, tol=1e-1, max_step=200)
n_arr = find_point(check_func, r0_guess, tol=1e-6)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, n_arr, label="二分数值求解")
ax.plot(r_arr,guess_func(r_arr), "--", label=r"猜测的解析解 $(\sqrt{4r+1}+1)/2$")
ax.grid()
ax.set_xlabel(r"$r$")
ax.set_ylabel(r"吸引域上界 $N$")
ax.set_title("稳定周期一轨道内吸引域上界随 $r$ 变化")
ax.legend()
display_inline(fig)

可以发现在整个周期一范围内数值求解的吸引域和我们估计的解析解基本吻合.

周期二

解析求解

类似周期一，定义 f2：

sym.var("r,x") # 定义符号
f2 = fn(r,x,2)
display_math(r"f^2(x)=",f2)

$\displaystyle f^2(x)=r - \left(r - x^{2}\right)^{2}$

计算 $f^2(x)$ 的不动点 $x^*$，并计算 $(f^2)'(x)$，将 $x^*$ 带入 $(f^2)'(x)$：

sol = sym.solve(f2-x,x)
df2 = sym.diff(f2,x)

display_math(r"\text{函数 } f^2(x)=",f2.simplify(),f"\\text{{ 一共找到 {len(sol)} 个不动点：}}")
args_list = []
for s in sol:
    args_list.extend([r"x^*=", s, r",\quad "])
display_math(*args_list)

display_math(r"f^2(x) \text{ 的导数为 }(f^2)'(x)=",df2,r"\text{, 带入 } x^* \text{ 有：}")
for s in sol:
    display_math(r"(f^2)'\left(", s, r"\right)=", df2.subs(x,s).simplify())

$\displaystyle \text{函数 } f^2(x)=r - \left(r - x^{2}\right)^{2}\text{ 一共找到 4 个不动点：}$

$\displaystyle x^*=\frac{1}{2} - \frac{\sqrt{4 r - 3}}{2},\quad x^*=\frac{\sqrt{4 r - 3}}{2} + \frac{1}{2},\quad x^*=- \frac{\sqrt{4 r + 1}}{2} - \frac{1}{2},\quad x^*=\frac{\sqrt{4 r + 1}}{2} - \frac{1}{2},\quad $

$\displaystyle f^2(x) \text{ 的导数为 }(f^2)'(x)=4 x \left(r - x^{2}\right)\text{, 带入 } x^* \text{ 有：}$

$\displaystyle (f^2)'\left(\frac{1}{2} - \frac{\sqrt{4 r - 3}}{2}\right)=4 - 4 r$

$\displaystyle (f^2)'\left(\frac{\sqrt{4 r - 3}}{2} + \frac{1}{2}\right)=4 - 4 r$

$\displaystyle (f^2)'\left(- \frac{\sqrt{4 r + 1}}{2} - \frac{1}{2}\right)=4 r + 2 \sqrt{4 r + 1} + 2$

$\displaystyle (f^2)'\left(\frac{\sqrt{4 r + 1}}{2} - \frac{1}{2}\right)=4 r - 2 \sqrt{4 r + 1} + 2$

$f^2(x)$ 的不动点一共有四个，可以分成下面两组：

$$ x^* = \frac{1\pm\sqrt{4r-3}}{2}, \quad x^* = \frac{\pm\sqrt{4r+1}-1}{2} $$

其中第二组就是之前我们求解出的周期一轨道，显然周期一轨道也是周期二轨道. 而第一组是新出现的周期二轨道. 且第一组的两个不动点的稳定条件相同，都是 $|4-4r|<1$，下面我们分别讨论每一个不动点的稳定性.

要使得不动点： $$ x^* = \frac{1\pm\sqrt{4r-3}}{2} $$ 稳定，则 $|4-4r|<1$，解得：$\frac{3}{4} < r < \frac{5}{4}$. 因此，当 $\frac{3}{4} < r < \frac{5}{4}$ 时，系统有稳定的周期二轨道：

$$ \dots \to \left( \frac{1+\sqrt{4r-3}}{2} \to \frac{1-\sqrt{4r-3}}{2} \right) \to \left( \frac{1+\sqrt{4r-3}}{2} \to \frac{1-\sqrt{4r-3}}{2} \right) \to \dots $$

---

要使得不动点：

$$ x^* = \frac{-\sqrt{4r+1}-1}{2} $$

稳定，则：

$$ |4r + 2\sqrt{4r+1}+2|<1 $$

由于 $r>0$，显然上式无解，因此此不动点不稳定.

---

要使得不动点：

$$ x^* = \frac{\sqrt{4r+1}-1}{2} $$

稳定，则：

$$ |4r - 2\sqrt{4r+1}+2|<1 $$

解得

sym.reduce_inequalities([df2.subs(x,sol[3])>-1, df2.subs(x,sol[3])<1], r)

$\displaystyle - \frac{1}{4} \leq r \wedge - \frac{1}{4} < r \wedge r < \frac{3}{4} \wedge r < \infty$

也就是 $0<r<\frac{3}{4}$（注意 $r>0$），和周期一轨道的稳定性相同（显然稳定的周期一轨道就是稳定的周期二轨道），并没有显现出和周期一不同的性质，因此这里不计入周期二轨道的范围内.

---

综上，当 $\frac{3}{4} < r < \frac{5}{4}$ 时，系统有稳定的周期二轨道：

$$ \dots \to \left( \frac{1+\sqrt{4r-3}}{2} \to \frac{1-\sqrt{4r-3}}{2} \right) \to \left( \frac{1+\sqrt{4r-3}}{2} \to \frac{1-\sqrt{4r-3}}{2} \right) \to \dots $$

数值验证

给定 $r=1$，我们在 $x\in[-2,2]$ 之间均匀采样 $5$ 个点，绘制其迭代过程：

fig, ax = plt.subplots(figsize=(8,6))

n = 20
r = 1
with np.errstate(over='ignore', invalid='ignore'):
    ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(-2,2,5), r, n, markersize=4, marker='o', linestyle='--')
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,21,2))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n})")
ax.axhline((1+np.sqrt(4*r-3))/2, color='c', linestyle=':', label=r"预期值1：$(1+\sqrt{4r-3})/2$")
ax.axhline((1-np.sqrt(4*r-3))/2, color='m', linestyle=':', label=r"预期值2：$(1-\sqrt{4r-3})/2$")
ax.grid()
ax.legend()

display_inline(fig)

可以发现同样有部分范围的 $x_0$ 发散到 $-\infty$，我们选取更多的 $x_0$ 并将吸引到周期二轨道的轨迹染成绿色：

fig, ax = plt.subplots(figsize=(8,6))

n = 50
r = 1.0
x0 = np.linspace(-2,2,50)
map_func = lambda r,x: r-x**2
with np.errstate(over='ignore', invalid='ignore'):
    n_arr = np.arange(n+1)
    x_arr = np.zeros((n+1, len(x0)), dtype=np.float64)
    x_arr[0] = x0.copy()
    for i in range(n):
        x_arr[i+1] = map_func(r, x_arr[i])
    for i, x in enumerate(x0):
        line_color = "red" if np.isinf(x_arr[-1,i]) else "green"
        ax.plot(n_arr, x_arr[:,i], label=f"$x_0$={x:.3f}",  markersize=4, marker='o', linestyle='--', color = line_color)
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,n+1,5))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，50 个初始值，绿线被吸引)")
ax.axhline((1+np.sqrt(4*r-3))/2, color='c', linestyle=':', label=r"预期值1：$(1+\sqrt{4r-3})/2$")
ax.text(6, (1+np.sqrt(4*r-3))/2 + 0.1, r"稳定的周期二轨道值1：$(1+\sqrt{4r-3})/2\approx$"+f"{(1+np.sqrt(4*r-3))/2 : .3f}", color='c', size=12)
ax.axhline((1-np.sqrt(4*r-3))/2, color='m', linestyle=':', label=r"预期值2：$(1-\sqrt{4r-3})/2$")
ax.text(6, (1-np.sqrt(4*r-3))/2 - 0.2, r"稳定的周期二轨道值2：$(1-\sqrt{4r-3})/2\approx$"+f"{(1-np.sqrt(4*r-3))/2 : .3f}", color='m', size=12)
ax.axhline(1.65, color='blue', linestyle=':')
ax.text(12, 1.65+0.1, r"分界线 1, $x_0\approx=1.65$", color='b', size=12)
ax.axhline(-1.65, color='blue', linestyle=':')
ax.text(12, -1.65-0.2, r"分界线 2, $x_0\approx-1.65$", color='b', size=12)
ax.grid()
#ax.legend()

display_inline(fig)

同样可以看到，类似周期一的情况，当初始值在两条分界线之间时，$x_n$ 收敛到周期二轨道，反之到达 $-\infty$.

对于周期二，我们除了可以关注是否被周期二轨道吸引外，还可以关注是吸引到 $p_1 \to p_2 \to \dots $ 还是 $p_2 \to p_1 \to \dots $，其中 $p_{1,2}$ 是周期二轨道的两个值：

$$ p_{1,2} = \frac{1+\pm\sqrt{4r-3}}{2} $$

我们根据最后一次迭代是靠近 $p_1$ 还是靠近 $p_2$ 来染色：

fig, ax = plt.subplots(figsize=(8,6))

n = 50
r = 1.0
x0 = np.linspace(-2,2,50)
map_func = lambda r,x: r-x**2
x_star1 = (1 + np.sqrt(4*r-3)) / 2.0
x_star2 = (1 - np.sqrt(4*r-3)) / 2.0
with np.errstate(over='ignore', invalid='ignore'):
    n_arr = np.arange(n+1)
    x_arr = np.zeros((n+1, len(x0)), dtype=np.float64)
    x_arr[0] = x0.copy()
    for i in range(n):
        x_arr[i+1] = map_func(r, x_arr[i])
    for i, x in enumerate(x0):
        x_f = x_arr[-1,i]
        line_color = "red" if np.isinf(x_arr[-1,i]) else "green"
        if np.isinf(x_f):
            line_color = "red"
        else:
            line_color = "c" if (np.abs(x_f - x_star1) - np.abs(x_f - x_star2) < 0) else "m"
        ax.plot(n_arr, x_arr[:,i], label=f"$x_0$={x:.3f}",  markersize=4, marker='o', linestyle='--', color = line_color)
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,n+1,5))
ax.set_ylabel("$x_n$")
ax.set_ylim((-2.5,2.5))
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，50 个初始值)")
ax.axhline((1+np.sqrt(4*r-3))/2, color='c', linestyle=':', label=r"预期值1：$(1+\sqrt{4r-3})/2$")
ax.text(6, (1+np.sqrt(4*r-3))/2 + 0.1, r"稳定的周期二轨道值1：$(1+\sqrt{4r-3})/2\approx$"+f"{(1+np.sqrt(4*r-3))/2 : .3f}", color='c', size=12)
ax.axhline((1-np.sqrt(4*r-3))/2, color='m', linestyle=':', label=r"预期值2：$(1-\sqrt{4r-3})/2$")
ax.text(6, (1-np.sqrt(4*r-3))/2 - 0.2, r"稳定的周期二轨道值2：$(1-\sqrt{4r-3})/2\approx$"+f"{(1-np.sqrt(4*r-3))/2 : .3f}", color='m', size=12)
ax.axhline(1.65, color='blue', linestyle=':')
ax.text(12, 1.65+0.1, r"分界线 1, $x_0\approx=1.65$", color='b', size=12)
ax.axhline(-1.65, color='blue', linestyle=':')
ax.text(12, -1.65-0.2, r"分界线 2, $x_0\approx-1.65$", color='b', size=12)
ax.grid()
#ax.legend()

display_inline(fig)

可以看到，在整体的周期二的吸引域内，不同轨道的吸引域呈现某种更精细的结构.

下面我们尝试求解周期二吸引域的范围. 并计算不同轨道的吸引域.

解析估计吸引域

对于整个周期二轨道的吸引域，类似前面的周期一，我们知道，若 $x_0 \in [-M,M]$ 其中：

$$ M = \frac{\sqrt{4r+1}+1}{2} $$

则在整个迭代范围中 $x_n \in [-M,M]$，而前面的数值计算中我们发现对于不同的 $x_0$，其要么趋向 $-\infty$，要么吸引到周期二轨道，因此，我们可以猜测整个周期二轨道的吸引域也为 $[-M,M]$

（如果周期一的吸引域是 $[-M,M]$，那么周期二的吸引域也应为 $[-M,M]$，因为在证明（虽然我没有能力证明）周期一的吸引域是 $[-M,M]$ 等价于证明 $f^n(x)$ 在这个范围内是压缩映射，也因此 $\left(f^{2}\right)^{n/2} (x)$ 也是压缩映射）.

而对于周期二吸引域内两条轨道单独的吸引域，这里暂时没有什么想法进行解析计算或者估计. 如果之后后想法再在这里补充.

解析求解过于困难，下面我们尝试通过数值计算寻找吸引域（尤其是两条轨道各自的吸引域）

数值计算吸引域

首先我们先验证整个周期二轨道的吸引域是否和我们估计的相同：

def make_check_func(r_arr, tol=1e-3, max_step=100): # 生成判断函数，向量化
    r_arr = np.array(r_arr,ndmin=1)
    x_star1 = (1 + np.sqrt(4*r_arr-3)) / 2.0
    x_star2 = (1 - np.sqrt(4*r_arr-3)) / 2.0
    def check_func(x0):
        x = x0.copy()
        with np.errstate(over='ignore', invalid='ignore'):
            for _ in range(max_step):
                x = r_arr - x**2
        return (np.abs(x - x_star1) < tol) | (np.abs(x - x_star2) < tol) 
    return check_func

n_r = 1000
guess_func = lambda r: (np.sqrt(4*r+1)+1)/2
r_arr = np.linspace(3/4,5/4,n_r)
r0_guess = np.array([[0,2] * n_r])
check_func = make_check_func(r_arr, tol=1e-1, max_step=300)
n_arr = find_point(check_func, r0_guess, tol=1e-6)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, n_arr, label="二分数值求解")
ax.plot(r_arr,guess_func(r_arr), "--", label=r"猜测的解析解 $(\sqrt{4r+1}+1)/2$")
ax.grid()
ax.set_xlabel(r"$r$")
ax.set_ylabel(r"吸引域上界 $N$")
ax.set_title("稳定周期二轨道内吸引域上界随 $r$ 变化")
ax.legend()
display_inline(fig)

可以看到，在整个稳定的周期二轨道的范围内，数值计算的周期二轨道的吸引域和我们估计的相同，即 $[-M,M]$，且：

$$ M = \frac{\sqrt{4r+1}+1}{2} $$

下面我们寻找整个周期二吸引域内两条轨道各自的吸引域. 由于对称性，显然吸引域是关于 $0$ 对称的，因此下面我们只计算大于 $0$ 的分界点.

和整个周期二吸引域只有一个分界点不同，两条轨道的分界点会有多个，因此这里我们采用如下方法寻找分界点.

这里我们对 $x_0 \in [0,M]$ 进行等间距采样（上界到 $M$ 保证一定会收敛到周期二而不是发散），并迭代一定次数，并判断最后一次迭代收敛到哪一个点，到 $p_1$ 记为 $1$，到 $p_2$ 记为 $0$，从而得到一个数字串：$111122222\dots$，找到里面每一个 $1,2$ 分界处，分界处的两个值的平均就是分界点的估计.

如果需要进一步细化也可以用每个分界处的两个值作为二分查找的初始区间，这里暂时先不用二分查找细化.

@nb.njit()
def classify(r_arr, n_grid=100, max_step=500):
    """
    对所有参数 r 和吸引域初值做迭代，返回分类结果
    """
    x_star1 = (1 + np.sqrt(4*r_arr-3))/2
    x_star2 = (1 - np.sqrt(4*r_arr-3))/2

    nr = len(r_arr)

    x_min = 0
    x_max = (np.sqrt(4*r_arr+1)+1)/2
    xs = np.linspace(-1, 1, n_grid)
    x0 = (x_max-x_min)[:,None]/2 * xs[None,:] + (x_max+x_min)[:,None]/2
    x = x0.copy()
    for _ in range(max_step):
        x = r_arr[:,None] - x**2

    d1 = np.abs(x - x_star1[:,None])
    d2 = np.abs(x - x_star2[:,None])
    labels = np.where(d2<d1, 1, 0)
    return x0, labels

def find_boundaries(r_arr, n_grid=100, max_step=500):
    x0s, labels = classify(r_arr, n_grid=n_grid, max_step=max_step)
    result = []
    for x0, label, r in zip(x0s, labels, r_arr):
        idxs = np.where(np.diff(label) != 0)[0]
        r_bnd = []
        for idx in idxs:
            r_bnd.append((x0[idx]+x0[idx+1])/2)
        result.append(r_bnd)
    return result

r_arr = np.linspace(3/4+1e-5, 5/4-1e-5, 1000)
sol = find_boundaries(r_arr, n_grid=500)
counts_sol = np.array([len(s) for s in sol])
total = counts_sol.sum()
x_all = np.repeat(r_arr, counts_sol)
y_all = np.concatenate(sol)

fig, ax = plt.subplots(figsize=(8,6))
bnd_func = lambda r : (np.sqrt(4*r+1)+1)/2
ax.plot(r_arr, bnd_func(r_arr), 'r', label="整个周期二吸引域分界")
ax.scatter(x_all, y_all, s=1, c="blue", label="不同轨迹吸引域的分界")
ax.grid()
ax.set_xlabel(r"$r$")
ax.set_ylabel(r"吸引域分界 $x_0$")
ax.set_title("稳定周期二轨道内吸引域分界随 $r$ 变化")
ax.legend()
display_inline(fig)

可以发现，在靠近整个周期二吸引域的分界处，不同轨道的吸引域变化很快，我们在附近放大绘制：

@nb.njit()
def classify2(r_arr, n_grid=100, max_step=500):
    """
    对所有参数 r 和吸引域初值做迭代，返回分类结果
    """
    x_star1 = (1 + np.sqrt(4*r_arr-3))/2
    x_star2 = (1 - np.sqrt(4*r_arr-3))/2

    nr = len(r_arr)

    x_min = 0.95 * (np.sqrt(4*r_arr+1)+1)/2
    x_max = (np.sqrt(4*r_arr+1)+1)/2
    xs = np.linspace(-1, 1, n_grid)
    x0 = (x_max-x_min)[:,None]/2 * xs[None,:] + (x_max+x_min)[:,None]/2
    x = x0.copy()
    for _ in range(max_step):
        x = r_arr[:,None] - x**2

    d1 = np.abs(x - x_star1[:,None])
    d2 = np.abs(x - x_star2[:,None])
    labels = np.where(d2<d1, 1, 0)
    return x0, labels

def find_boundaries(r_arr, n_grid=100, max_step=500):
    x0s, labels = classify2(r_arr, n_grid=n_grid, max_step=max_step)
    result = []
    for x0, label, r in zip(x0s, labels, r_arr):
        idxs = np.where(np.diff(label) != 0)[0]
        r_bnd = []
        for idx in idxs:
            r_bnd.append((x0[idx]+x0[idx+1])/2)
        result.append(r_bnd)
    return result

r_arr = np.linspace(3/4+1e-5, 5/4-1e-5, 1000)
sol = find_boundaries(r_arr, n_grid=500)
counts_sol = np.array([len(s) for s in sol])
total = counts_sol.sum()
x_all = np.repeat(r_arr, counts_sol)
y_all = np.concatenate(sol)
fig, ax = plt.subplots(figsize=(8,6))
bnd_func = lambda r : (np.sqrt(4*r+1)+1)/2
ax.plot(r_arr, bnd_func(r_arr), 'r', label="整个周期二吸引域分界")
ax.scatter(x_all, y_all, s=1, c="blue", label="不同轨迹吸引域的分界")
ax.grid()
ax.set_xlabel(r"$r$")
ax.set_ylabel(r"吸引域分界 $x_0$")
ax.set_title("稳定周期二轨道内吸引域分界随 $r$ 变化")
ax.legend()
display_inline(fig)

可以看到，越接近这个周期二吸引域的边界，不同轨道的吸引域的分界越密集.

我们可以给图像绘制颜色显示不同初值的最后的去向：

n = 100
x0_arr = np.linspace(0, 2, 1000)
r_arr = np.linspace(3/4+1e-5, 5/4-1e-5, 1000)
X0, Y0 = np.meshgrid(r_arr, x0_arr)
with np.errstate(over='ignore', invalid='ignore'):
    x_star1 = (1 + np.sqrt(4*X0-3)) / 2.0
    x_star2 = (1 - np.sqrt(4*X0-3)) / 2.0
    Y = Y0.copy()
    for _ in range(n):
        Y = X0 - Y**2
result = np.zeros(X0.shape, dtype=int)
result[np.isnan(Y)] = 0
result[~np.isnan(Y) & (np.abs(Y-x_star1) < np.abs(Y-x_star2))] = 1
result[~np.isnan(Y) & (np.abs(Y-x_star1) > np.abs(Y-x_star2))] = 2

colors = ['red', 'c', 'm']
cmap = mpl.colors.ListedColormap(colors)

fig, ax = plt.subplots(figsize=(8,6))
ax.imshow(result, cmap=cmap, origin='lower',
           extent=[r_arr.min(), r_arr.max(), x0_arr.min(), x0_arr.max()],
           aspect='auto')
ax.set_xlabel("$r$")
ax.set_ylabel("$x_0$")
ax.set_title(r"不同条件下的渐进行为（红：到负无穷，青：$p_1 \to p_2$ 循环，洋红：$p_2 \to p_1$ 循环）")
display_inline(fig)

在 $1.5$ 附近放大：

n = 100
x0_arr = np.linspace(1.5-0.1, 1.5+0.1, 1000)
r_arr = np.linspace(3/4+1e-5, 5/4-1e-5, 1000)
X0, Y0 = np.meshgrid(r_arr, x0_arr)
with np.errstate(over='ignore', invalid='ignore'):
    x_star1 = (1 + np.sqrt(4*X0-3)) / 2.0
    x_star2 = (1 - np.sqrt(4*X0-3)) / 2.0
    Y = Y0.copy()
    for _ in range(n):
        Y = X0 - Y**2
result = np.zeros(X0.shape, dtype=int)
result[np.isnan(Y)] = 0
result[~np.isnan(Y) & (np.abs(Y-x_star1) < np.abs(Y-x_star2))] = 1
result[~np.isnan(Y) & (np.abs(Y-x_star1) > np.abs(Y-x_star2))] = 2
colors = ['red', 'c', 'm']
cmap = mpl.colors.ListedColormap(colors)

fig, ax = plt.subplots(figsize=(8,6))
ax.imshow(result, cmap=cmap, origin='lower',
           extent=[r_arr.min(), r_arr.max(), x0_arr.min(), x0_arr.max()],
           aspect='auto')
ax.set_xlabel("$r$")
ax.set_ylabel("$x_0$")
ax.set_title(r"不同条件下的渐进行为（红：到负无穷，青：$p_1 \to p_2$ 循环，洋红：$p_2 \to p_1$ 循环）")
display_inline(fig)

综上，我们知道了周期二轨道的整体的吸引域为 $[-M,M]$，其中：

$$ M = \frac{\sqrt{4r+1}+1}{2} $$

且在整体的吸引域内各条轨迹的吸引域具有丰富细节，且越靠近整体吸引域的边界变化越快.

为何会出现 Overflow

在前面的绘制分叉图的过程中，我们遇到了 Overflow 的问题，具体而言，Overflow 发生在当数值非常大或者非常小以至于超过计算机中的表达范围时（在 Python 中默认使用 float64，其最大能保存约 $10^{-308}$ 到 $10^{308}$ 范围内的数值（精度约 16 位）），超过范围的数会变成 inf 或者 -inf 从而无法进一步进行数值计算，并出现 Overflow 问题.

回顾我们绘制分叉图的过程，我们多次迭代计算了 $f(x)=r-x^2$，这里 $r$ 是一个比较小的值（$r\in[0,2]$），因此出现 Overflow 只能出现在计算 $x^2$ 之后. 我们来尝试分析给定定义域 $x\in[-2,2]$ 经过 $f^n (x)$ 后值域 $\left\{ f^n(x) | x \in [-2,2] \right\}$ 中的元素绝对值的最大值（也就是 $\text{max} \left\{ |f^n(x)| \right\} $）的变化：

由于 $f(x)$ 是关于 $x=0$ 对称的开口向下的抛物线，因此 $\text{max} \left\{ |f(x)| \right\} $ 要么出现在 $f(0)$ 要么出现在 $f(2)$（由对称性 $f(-2)=f(2)$，因此不再考虑 $-2$，下同），由于 $\left\{ f(x) \right\}$ 一定是 $[-\text{max} \left\{ |f(x)| \right\},\text{max} \left\{ |f(x)| \right\}]$ 的子集，即：

$$ \left\{ f(x) \right\} \subset [-\text{max} \left\{ |f(x)| \right\},\text{max} \left\{ |f(x)| \right\}] $$

因此 $\text{max} \left\{ |f^2(x)| \right\} $ 要么出现在 $|f(0)|$，要么出现在 $\left| f\left( \text{max} \left\{ |f(x)| \right\} \right) \right|$，后者等价于 $\text{max} \left\{ \left| f^2(0)\right|, \left|f^2(2) \right| \right\}$.

上述过程可以推到 $f^n$，综上：

$$ \text{max} \left\{ |f^n(x)| \right\} = \text{max} \left\{ |f(0)|, \left| f^n (0) \right|, \left| f^n (2) \right| \right\}$$

而由于 $|f^n(2)|$ 当 $r<2$ 时随着 $n$ 增大而增大（比较显然，不写详细证明了），直至 $\infty$，因此，我们知道 $\text{max} \left\{ |f^n(x)| \right\} $ 一定会随着 $n$ 增大而增大，直至到达双精度浮点表示的最大值，从而出现 Overflow.

全周期分叉图

分叉图

绘制全周期分叉图的方法和第三次作业类似，这里不再进行说明，直接使用上次作业的函数（进行了一点优化）：

def bifurcation_diagram(ax, map_func, r_range, x_range, n_warm, n_sample, **kwargs):
    r = np.linspace(r_range[0], r_range[1], r_range[2])
    x = np.random.uniform(x_range[0], x_range[1], x_range[2])
    r_tmp = r.reshape(r_range[2], 1) # 转为列向量，为了利用 broadcast

    # 到达稳定状态
    for _ in range(n_warm):
        x = map_func(r_tmp, x)

    result = np.zeros((r_range[2], x_range[2], n_sample))
    for i in range(n_sample):
        x = map_func(r_tmp, x)
        result[:, :, i] = x
    x_plot = result.reshape(r_range[2], -1)
    ax.plot(r, x_plot, **kwargs)
    ax.set_ylim((x_range[0], x_range[1]))
    ax.set_xlim((r_range[0], r_range[1]))
    return ax, r, x_plot

def bifurcation_diagram_from_save(ax, r_save, x_save, r_lim, x_lim, **kwargs):
    ax.plot(r_save, x_save, **kwargs)
    ax.set_ylim((x_lim[0], x_lim[1]))
    ax.set_xlim((r_lim[0], r_lim[1]))
    return ax

在 $r\in[0,2]$ 内均匀采样 $5 \times 10^3$ 个点，$x\in[-2,2]$ 内均匀随机选取 $20$ 个点，记录最后 $100$ 次的值（一个 $r$ 一共产生 $2\times 10^3$ 个点）.

这里我们忽略最终到达 $\pm \infty$ 的点（也就是我们只记录被吸引的点）

fig, ax = plt.subplots(figsize=(8,6))
with np.errstate(over='ignore', invalid='ignore'):
    ax, r_save, x_save = bifurcation_diagram(ax, lambda r,x: r-x**2, (0,2,5000), (-2,2,20), 1000, 10, marker=',', color='k', linestyle='none', alpha=0.1)
ax.set_xlabel("$r$")
ax.set_ylabel("$x_n$")
ax.set_title(r"分叉图（$x_{n+1} = r-x^2$）")
ax.grid()
display_inline(fig)

可以看见我们得到了一个类似前一次作业的 Logistic Map 的分叉图.

周期三

数值确定周期三位置

从全周期分叉图中我们可以看到，大概在 $r=1.75$ 附近出现周期三轨道. 这里我们针对这附近的范围进行放大：

fig, ax = plt.subplots(figsize=(8,6))
with np.errstate(over='ignore', invalid='ignore'):
    ax, _, _ = bifurcation_diagram(ax, lambda r,x: r-x**2, (1.74,1.78,1000), (-2,2,20), 1000, 10, marker=',', color='k', linestyle='none', alpha=0.1)
ax.set_xlabel("$r$")
ax.set_ylabel("$x_n$")
ax.set_title(r"分叉图（$x_{n+1} = r-x^2$，$r\in[1.74,1.78]$）")
ax.axvline(1.75, color='blue', linestyle=':')
ax.text(1.7502, -1.75, r"$r=1.75$", color='b', size=12)
ax.grid()
display_inline(fig)

可以看到，确实在 $r=1.75$ 附近出现周期三轨道. 我们绘制迭代过程验证：

fig, ax = plt.subplots(figsize=(8,6))

n = 100
r = 1.75
ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(0,0,1), r, n, markersize=4, marker='o', linestyle='--')
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,n+1,10))
ax.set_ylabel("$x_n$")
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，$x_0$=0)")
ax.grid()
#ax.legend()

display_inline(fig)

解析求解周期三位置

实际上，周期三轨道的 $r$ 是可以通过解析求解的，我们绘制 $y=f^3(x)$ 和 $y=x$（取 $r=1.75$）

fig, ax = plt.subplots(figsize=(8,6))
x_arr = np.linspace(-2,2,300)
ax.plot(x_arr, fn(1.75, x_arr, 3), label=r"$y=f^3(x)$")
ax.plot(x_arr, x_arr, label="$y=x$")
ax.set_xlabel("$x$")
ax.set_ylabel("$y$")
ax.set_title("$y=f^3(x)$ 和 $y=x$")
ax.set_xlim((-2,2))
ax.set_ylim((-2,2))
ax.grid()
ax.legend()
display_inline(fig)

可以发现，要出现周期三，需要 $f^3(x)$ 和 $x$ 出现“相切”，因此，要出现周期三，需要满足以下两个方程：

• $x^* = f^3(x^*)$
• $(f^3)'(x^*)=1$

两个方程解两个未知数（$x^*$ 和 $r$），因此是可以解的，下面使用 sympy 给出解析解：

定义两个方程：

sym.var("x,r")
f3 = fn(r,x,3)
df3 = sym.diff(f3,x)
P = f3 - x
Q = df3 - 1
display_math(P,r"=0")
display_math(Q,r"=0")

$\displaystyle r - x - \left(r - \left(r - x^{2}\right)^{2}\right)^{2}=0$

$\displaystyle - 8 x \left(r - x^{2}\right) \left(r - \left(r - x^{2}\right)^{2}\right) - 1=0$

对于多项式，根据数学上有关多项式的理论，我们可以通过计算结式（Resultant）来消去 $x$，从而得到关于 $r$ 的多项式方程：

tmp = sym.resultant(P,Q,x)
display_math(r"0=",tmp.expand())

$\displaystyle 0=- 16777216 r^{12} + 100663296 r^{11} - 251658240 r^{10} + 384827392 r^{9} - 449839104 r^{8} + 402653184 r^{7} - 275767296 r^{6} + 160161792 r^{5} - 59006976 r^{4} + 21512960 r^{3} + 823543$

因式分解，有：

display_math(r"0=",sym.factor(tmp))

$\displaystyle 0=- \left(4 r - 7\right)^{3} \left(4 r + 1\right) \left(16 r^{2} - 4 r + 7\right)^{4}$

显然，上面三个因式有且只有第一个有大于 $0$ 的实根，根为：

$$ r = \frac{7}{4} = 1.75 $$

由此我们证明了当出现周期三时，$r=7/4$.

注：上述过程对于 Logistic Map 同样适用，可以解得对于原始的 Logistic Map：$x_{n+1} = r x_n (1-x_n) $，其周期三轨道出现在 $r=1+2\sqrt{2}$

周期五

从全周期分叉图中我们可以看到，大概在 $r=1.625$ 附近出现周期五轨道. 这里我们针对这附近的范围进行放大：

fig, ax = plt.subplots(figsize=(8,6))
with np.errstate(over='ignore', invalid='ignore'):
    ax, _, _ = bifurcation_diagram(ax, lambda r,x: r-x**2, (1.615,1.635,1000), (-2,2,20), 1000, 10, marker=',', color='k', linestyle='none', alpha=0.1)
ax.set_xlabel("$r$")
ax.set_ylabel("$x_n$")
ax.set_title(r"分叉图（$x_{n+1} = r-x^2$，$r\in[1.615,1.635]$）")
ax.axvline(1.625, color='blue', linestyle=':')
ax.text(1.6252, -1.75, r"$r=1.625$", color='b', size=12)
ax.grid()
display_inline(fig)

可以看到，确实在 $r=1.625$ 附近出现周期五轨道. 我们绘制迭代过程验证：

fig, ax = plt.subplots(figsize=(8,6))

n = 100
r = 1.625
ax = plot_orbit(ax, lambda r,x:r-x**2, np.linspace(0,0,1), r, n, markersize=4, marker='o', linestyle='--')
ax.set_xlabel("$n$")
ax.set_xticks(np.arange(0,n+1,10))
ax.set_ylabel("$x_n$")
ax.set_title(f"$x_n$ 随 $n$ 变化过程 ($r$={r:.3f}, $n$={n}，$x_0$=0)")
ax.grid()
#ax.legend()

display_inline(fig)

可以发现，当 $r=1.625$ 时确实出现了周期五轨道.

全参数范围的李雅普诺夫指数

计算思路

根据李雅普诺夫指数的定义，对于微扰 $\delta_0$，迭代 $n$ 次后，有：

$$ \delta_n = \delta_0 \mathrm{e}^{\lambda n} $$

因此：

$$ \lambda = \frac{1}{n} \ln \frac{\delta_n}{\delta_0} = \frac{1}{n} \ln |(f^n)'(x)| = \frac{1}{n} \sum_{k=1}^n |f'(x_k)|$$

从而得到计算思路，对于一个 $r$，经过 $n_{\text{warm}}$ 去除暂态后，迭代 $n$ 次，记录 $\ln|f'(x_k)|$ 的和 $S$，最后 $\lambda = S/n$

绘制

考虑到之后可能会对不同的迭代过程绘制全参数范围的李雅普诺夫指数，这里将上述的通过解析表达式的计算过程抽象成一个函数（为了加快速度，要求 f 和 df 是向量化的）.

def calc_lyapunov_1d(f, df, x0, r_arr, n_warm, n_sample, offset=1e-300):
    x0 = np.full_like(r_arr, x0)
    for _ in range(n_warm):
        x0 = f(r_arr, x0)
    sum_ln = np.zeros_like(r_arr)
    for _ in range(n_sample):
        x0 = f(r_arr, x0)
        sum_ln += np.log(np.abs(df(r_arr, x0)+offset)) # offset, 防止 log(0)
    return sum_ln / n_sample

在 $r\in[0,2]$ 内均匀采样 $1 \times 10^4$ 个点，取 $x_0 = 0.1$，记录最后 $200$ 次的值来计算李雅普诺夫指数. 并绘制计算得到的李雅普诺夫指数和 $r$ 的关系：

r_arr = np.linspace(0,2,10000)
lambda_r = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            1000,
                            200)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, lambda_r, label=r"$\lambda$")
ax.set_ylim((-2,1))
ax.set_xlabel("$r$")
ax.set_ylabel(r"$\lambda$")
ax.set_title("全参数范围的李雅普诺夫指数")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

这里我们和前面的分叉图一同显示：

fig, axs = plt.subplots(2, 1, sharex=True, figsize=(8,6))
ax1, ax2 = axs

ax1 = bifurcation_diagram_from_save(ax1, r_save, x_save, (0,2), (-2,2), marker=',', color='k', linestyle='none', alpha=0.1)

ax1.set_xlabel("$r$")
ax1.set_ylabel("$x_n$")
ax1.set_title(r"分叉图和全参数范围的李雅普诺夫指数（$x_{n+1} = r-x^2$）")
ax1.grid()

ax2.plot(r_arr, lambda_r, label=r"$\lambda$")
ax2.grid()
ax2.set_ylim((-2,1))
ax2.set_xlabel("$r$")
ax2.set_ylabel(r"$\lambda$")
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
display_inline(fig)

可以发现，李雅普诺夫指数 $\lambda=0$ 时开始分叉，$\lambda>0$ 时出现混沌，当 $\lambda$ 很小时（实际上应该趋向于 $-\infty$，但是由于计算和绘图限制，没有显示）为超稳定轨道（对应 $x_n = 0$ 在最后的轨道中）.

倍周期分叉点、$r_{\infty}$、费根鲍姆常数

倍周期分叉点

根据前面的李雅普诺夫指数和分叉图的图可以发现，对于每一个倍周期分叉点，其李雅普诺夫指数都有 $\lambda=0$. 因此，我们可以搜索全参数范围的李雅普诺夫指数中约等于 $0$ 的点.

r_arr[np.abs(lambda_r) < 2e-3]

array([0.74827483, 0.74847485, 0.74867487, 0.74887489, 0.74907491,
       0.74927493, 0.74947495, 0.74967497, 0.74987499, 0.75007501,
       0.75027503, 0.75047505, 0.75067507, 1.24912491, 1.24932493,
       1.24952495, 1.24972497, 1.24992499, 1.25012501, 1.25032503,
       1.36773677, 1.36793679, 1.36813681, 1.39393939, 1.39413941,
       1.40114011, 1.40134013, 1.40194019, 1.52215222])

可以看到，受限于计算精度，在设置允差为 $2 \times 10^{-3}$ 时我们找到了一堆可能的点. 当然我们可以通过更精细地设置允差来给出更好的筛选，但这个过程太依赖于我们的人工观察和经验（判断允差是否合适使得没有多余点也没有少点是依赖于人的观察的）.

当然，上述结果可以尝试再通过聚类等方式进行优化. 这里则使用另外一种方式更好地找到倍周期分叉点：

初筛

首先，比较麻烦的是如何给定包含且只包含一个倍周期分叉点的小区间作为初始区间，这里我们采用如下办法：

1. 在某个范围内均匀取 $r$.
2. 对于某个 $r$，经过 $n_{warm}$ 后，开一个数组
3. 记录 $x_n$，如果数组中没有（距离数组中的元素小于 tol），则添加到数组中，并再迭代一次，否则认为出现循环，得到周期数
4. 对每个 $r$ 执行上面 2,3 两步，从而得到周期数，并给出周期数变化的位置作为初始值.

从全参数范围的李雅普诺夫指数中，我们可以知道，混沌区大约在 $r\approx1.4$ 附近取到，为了保险起见，我们放大在 $1.4$ 附近的李雅普诺夫指数：

r_min = 1.4010
r_max = 1.4012
r_arr = np.linspace(r_min,r_max,5000)
lambda_r = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            500,
                            10000,
                            offset=0)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, lambda_r, label=r"$\lambda$")
ax.set_xlabel("$r$")
ax.set_ylabel(r"$\lambda$")
ax.set_title(f"部分参数范围的李雅普诺夫指数（$r\\in[{r_min:.4f}, {r_max:.4f}]$）")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

可以发现，混沌区大约在 $r\approx1.40115$ 出现. 我们取 $r\in[0,1.40115]$. 并写出如下函数得到 $r\in [0,1.40115]$ 的周期数随 $r$ 的变化关系：

为了加速，这里使用 numba 进行 jit 加速以及并行计算加速，第一次运行可能需要有几秒的编译时间

@nb.njit()
def calc_cycle_single(r, x0=0.0, tol=1e-6, n_warm=1000, max_period=32):
    """
    检测给定 r 的周期数（最大为 max_period）
    超过 max_period 或无循环则返回 -1
    """
    x = x0
    for _ in range(n_warm):
        x = r - x * x

    xs = np.empty(max_period + 1)
    xs[0] = x
    n_stored = 1

    for n in range(1, max_period + 1):
        x = r - x * x
        # 检查是否与已有点接近
        for i in range(n_stored):
            if abs(x - xs[i]) < tol:
                return n - i
        xs[n_stored] = x
        n_stored += 1
    return -1

@nb.njit(parallel=True)
def calc_cycles(r_arr, x0=0.0, tol=1e-6, n_warm=1000, max_period=32):
    n = len(r_arr)
    periods = np.empty(n)
    for i in nb.prange(n):
        periods[i] = calc_cycle_single(r_arr[i], x0, tol, n_warm, max_period)
    return periods

r_arr = np.linspace(0, 1.401155, 2000)
ncycles_arr = calc_cycles(r_arr, max_period=128, tol=1e-4, n_warm=10000)
mask = (ncycles_arr > 0)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr[mask], ncycles_arr[mask],'.')
for n in range(7):
    #ax.text(0.2, 2**(n)+0.1, f"$n=2^{n}$", color='k', size=12)
    ax.axhline(2**(n), linestyle=':', color='k')
ax.grid()
ax.set_xlabel("$r$")
ax.set_ylabel("周期数 $n$")
ax.set_title("倍周期范围内周期数随 $r$ 变化（黑色虚线为 $2^n$）")
display_inline(fig)

对于 $r\in[1.3,1.40115]$ 附近的点进一步细化：

r_arr2 = np.linspace(1.39, 1.40116, 2000)
ncycles_arr2 = calc_cycles(r_arr2, max_period=128, tol=1e-4, n_warm=10000)
mask2 = (ncycles_arr2 > 0)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr2[mask2], ncycles_arr2[mask2],'.')
for n in range(7):
    #ax.text(0.2, 2**(n)+0.1, f"$n=2^{n}$", color='k', size=12)
    ax.axhline(2**(n), linestyle=':', color='k')
ax.grid()
ax.set_xlabel("$r$")
ax.set_ylabel("周期数 $n$")
ax.set_title("倍周期范围内周期数随 $r$ 变化（黑色虚线为 $2^n$）")
display_inline(fig)

r_arr3 = np.linspace(1.400, 1.40116, 2000)
ncycles_arr3 = calc_cycles(r_arr3, max_period=128, tol=1e-6, n_warm=10000)
mask3 = (ncycles_arr3 > 0)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr3[mask3], ncycles_arr3[mask3],'.')
for n in range(8):
    #ax.text(0.2, 2**(n)+0.1, f"$n=2^{n}$", color='k', size=12)
    ax.axhline(2**(n), linestyle=':', color='k')
ax.grid()
ax.set_xlabel("$r$")
ax.set_ylabel("周期数 $n$")
ax.set_title("倍周期范围内周期数随 $r$ 变化（黑色虚线为 $2^n$）")
display_inline(fig)

从中我们可以得到每个倍周期分叉点的低精度估计：

def find_cycle_change(r_arr, ncycles_arr):
    result = []
    for i in range(len(r_arr) - 1):
        p1 = ncycles_arr[i]
        p2 = ncycles_arr[i + 1]
        if p1 <= 0 or p2 <= 0 or p1 == p2:
            continue
        if (p2 == 2 * p1):
            result.append(((r_arr[i] + r_arr[i + 1]) / 2, int(p1), int(p2)))
    return result

guess_pt = find_cycle_change(r_arr, ncycles_arr)[:3] + find_cycle_change(r_arr2, ncycles_arr2)[:2] + find_cycle_change(r_arr3, ncycles_arr3)
for (r, p1, p2) in guess_pt:
    print(f"{p1:<2} -> {p2:<3} 倍周期分叉点估计: {r:.6f}")

1  -> 2   倍周期分叉点估计: 0.749642
2  -> 4   倍周期分叉点估计: 1.250105
4  -> 8   倍周期分叉点估计: 1.367861
8  -> 16  倍周期分叉点估计: 1.394039
16 -> 32  倍周期分叉点估计: 1.399644
32 -> 64  倍周期分叉点估计: 1.400823
64 -> 128 倍周期分叉点估计: 1.401084

细化（0.618 法）

由于后面的费根鲍姆常数需要较高精度的倍周期分叉点，因此，这里需要进一步细化分叉点的值.

首先，前面我们得到了倍周期分叉点的估计，记为 $p_k$，由于这里的倍周期分叉点是不够精确的，因此这里需要使用二分法或者 0.618 法进一步细化，这就需要给出包含且只包含每个分叉点的小区间，这里按照如下方法定义：

对于分叉点 $p_k$，构造小区间为：

$$ \left( p_k - \frac{1}{2} (p_k - p_{k-1}), p_k + \frac{1}{2} (p_{k+1} - p_{k}) \right) $$

（头尾两个设半径为 0.001）从而得到：

guess_r = np.array([pt[0] for pt in guess_pt])

guess_r_range = np.zeros((7, 2))

for i in range(1,6):
    guess_r_range[i] = np.array([guess_r[i] - 0.5 * (guess_r[i] - guess_r[i-1]),
                                 guess_r[i] + 0.5 * (guess_r[i+1] - guess_r[i])])
guess_r_range[0] = np.array([guess_r[0]-0.001,
                             guess_r[0]+0.001])
guess_r_range[6] = np.array([guess_r[6]-0.001,
                             guess_r[6]+0.001])
for (r, p1, p2), (rr1, rr2) in zip(guess_pt, guess_r_range):
    print(f"{p1:<2} -> {p2:<3} 倍周期分叉点估计范围: {rr1:.6f} 到 {rr2:.6f}")

1  -> 2   倍周期分叉点估计范围: 0.748642 到 0.750642
2  -> 4   倍周期分叉点估计范围: 0.999874 到 1.308983
4  -> 8   倍周期分叉点估计范围: 1.308983 到 1.380950
8  -> 16  倍周期分叉点估计范围: 1.380950 到 1.396842
16 -> 32  倍周期分叉点估计范围: 1.396842 到 1.400234
32 -> 64  倍周期分叉点估计范围: 1.400234 到 1.400954
64 -> 128 倍周期分叉点估计范围: 1.400084 到 1.402084

对于每一个区间，我们知道，当位于倍周期分叉点时，$\lambda=0$，而不在倍周期分叉点时，$\lambda<0$，因此，倍周期分叉点在小区间内满足 $\lambda$ 为最值.

因此，我们将在小区间内寻找倍周期分叉点的问题转化为一个优化问题，且这个函数在每个区间内是“单峰的”，因此我们可以使用梯度上升法（通过差分近似微分估计梯度），0.618 法等多种数值方法来找最值.

这里我们使用 0.618 法，具体过程如下：

1. 给定边界 $[a,b]$，计算其中的两个黄金分割位置 $c,d$，其中 $c = a + (b-a)/\phi $，$d = b - (b-a)/\phi $，$\phi$ 为黄金分割比
2. 计算点 $c$ 和点 $d$ 位置的李雅普诺夫指数，并比较大小，若 $\lambda_{c} > \lambda_{d}$，则跳转第三步，否则跳转第四步
3. 令 $b \to d$，$d \to c$ 回到第一步
4. 令 $a \to c$，$c \to d$ 回到第一步
5. 重复上述过程，直到最后满足 $b-a < \text{tol}$

上述 0.618 法细化过程如下：

@nb.njit(parallel=True)
def refine_by_0618(r_range_arr, precision=1e-8, max_iter=5000):
    def lyapunov(r):
        x = 0.0
        for _ in range(10000): # 缓解临界慢化问题，n_warm 设置较大
            x = r - x**2
        sum_ln = 0.0
        for _ in range(10000):
            x = r - x**2
            sum_ln += np.log(2 * np.abs(x))
        return sum_ln / 10000

    phi = (np.sqrt(5) - 1) / 2
    n_r = r_range_arr.shape[0]
    results = np.zeros(n_r, dtype=np.float64)
    for i in nb.prange(n_r):
        a = r_range_arr[i, 0]
        b = r_range_arr[i, 1]
        c = b - phi * (b - a)
        d = a + phi * (b - a)
        f_c = lyapunov(c)
        f_d = lyapunov(d)

        for _ in range(max_iter):
            if (b - a) < precision:
                break
            if f_c > f_d:
                b = d
                d = c
                c = b - phi * (b - a)
                f_d = f_c
                f_c = lyapunov(c)
            else:
                a = c
                c = d
                d = a + phi * (b - a)
                f_c = f_d
                f_d = lyapunov(d)
        results[i] = (a + b) / 2
    return results

res = refine_by_0618(guess_r_range)

for (_, p1, p2), r in zip(guess_pt, res):
    print(f"{p1:<2} -> {p2:<3} 倍周期分叉点: {r:.6f}")

1  -> 2   倍周期分叉点: 0.749965
2  -> 4   倍周期分叉点: 1.249983
4  -> 8   倍周期分叉点: 1.368091
8  -> 16  倍周期分叉点: 1.394043
16 -> 32  倍周期分叉点: 1.399630
32 -> 64  倍周期分叉点: 1.400828
64 -> 128 倍周期分叉点: 1.400828

最后一个 64 到 128 的倍周期分叉点显然有问题，可能是因为我们的初始范围定的不好（考虑到最后一个点的初始范围是人为定的），这里暂时略去.

同时，对比我们已经知道的 1 到 2 分叉点（0.75），2 到 4 分叉点（1.25），上面的数值计算还是较为准确的.

因此，我们得到了数值计算的倍周期分叉点，如下：

res_rn = res[:-1].copy()
for (_, p1, p2), r in zip(guess_pt, res_rn):
    print(f"{p1:<2} -> {p2:<2} 倍周期分叉点: {r:.6f}")

1  -> 2  倍周期分叉点: 0.749965
2  -> 4  倍周期分叉点: 1.249983
4  -> 8  倍周期分叉点: 1.368091
8  -> 16 倍周期分叉点: 1.394043
16 -> 32 倍周期分叉点: 1.399630
32 -> 64 倍周期分叉点: 1.400828

$r_{\infty}$

初筛

在前面全参数范围的李雅普诺夫指数的图中，我们知道混沌取大约在 $r\approx1.40115$ 出现，我们在这个数邻域放大李雅普诺夫指数：

r_min = 1.401154
r_max = 1.401156
r_arr = np.linspace(r_min,r_max,5000)
lambda_r = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            1000,
                            10000,
                            offset=0)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, lambda_r, label=r"$\lambda$")
ax.set_xlabel("$r$")
ax.set_ylabel(r"$\lambda$")
ax.set_title(f"部分参数范围的李雅普诺夫指数（$r\\in[{r_min:.6f}, {r_max:.6f}]$）")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

更一般的，我们可以查询 $\lambda$ 何时第一次出现大于 $0$. （由于临界慢化，上面的一些分叉点的 $\lambda$ 略微小于 $0$，因此这里查询 $\lambda$ 第一次大于 $0$ 的 $r$ 实际高估了 $r_{\infty}$）

r_arr[np.where(lambda_r > 0)[0][0]]

$\displaystyle 1.40115516223245$

因此，一个 $r_{\infty}$ 的粗略估计为 $r_{\infty} = 1.401155$

细化（二分法）

同样，我们需要细化上面的计算，这里由于 $\lambda$ 在 $r_{\infty}$ 两侧变号，比较适合二分法，因此这里使用二分法进行细化，初始邻域宽度定为：$10^{-6}$

def find_r_inf(r_arr, tol=1e-10, max_step=50): # 二分查找函数，向量化
    # check_func 1 改变下界, 0 改变上界
    def check_func(r):
        x = 0.0
        for _ in range(10000): # 缓解临界慢化问题，n_warm 设置较大
            x = r - x**2
        sum_ln = 0.0
        for _ in range(10000):
            x = r - x**2
            sum_ln += np.log(2 * np.abs(x))
        lambda_r = sum_ln / 10000
        return lambda_r < 0
        
    a, b = r_arr[0], r_arr[1]
    for _ in range(max_step):
        mid = 0.5 * (a + b)
        check = check_func(mid)
        check = np.asarray(check, dtype=bool)
        a = np.where(check, mid, a)
        b = np.where(check, b, mid)
        if np.all(np.abs(b - a) < tol):
            break
    return 0.5 * (a + b)

r_inf_guess = 1.401155
result = find_r_inf(np.array([r_inf_guess - 1e-6, r_inf_guess + 1e-6]))
res_r_inf = result
res_r_inf

$\displaystyle 1.40115521878052$

因此，我们得到了数值计算的 $r_{\infty} = 1.40115522\dots$

验证倍周期分叉点和 $r_{\infty}$

我们将前面计算的倍周期分叉点和 $r_{\infty}$ 在分叉图中标准以验证：

r_arr = np.linspace(0,2,10000)
lambda_r = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            1000,
                            200)
fig, axs = plt.subplots(2, 1, sharex=True, figsize=(8,6))
ax1, ax2 = axs
ax1 = bifurcation_diagram_from_save(ax1, r_save, x_save, (0,2), (-2,2), marker=',', color='k', linestyle='none', alpha=0.1)
ax1.set_xlabel("$r$")
ax1.set_ylabel("$x_n$")
ax1.set_title(r"分叉图和全参数范围的李雅普诺夫指数（红：分叉点，绿：$r_{\infty}$）")
ax1.grid()

ax2.plot(r_arr, lambda_r, label=r"$\lambda$")
ax2.grid()
ax2.set_ylim((-2,1))
ax2.set_xlabel("$r$")
ax2.set_ylabel(r"$\lambda$")
ax2.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")

for ax in axs:
    for rn in res_rn:
        ax.axvline(rn, linestyle=":", color="r")
    ax.axvline(res_r_inf, linestyle=":", color="g")

display_inline(fig)

在接近混沌区放大：

r_arr = np.linspace(1.38,1.41,10000)
lambda_r = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            1000,
                            200)
fig, axs = plt.subplots(2, 1, sharex=True, figsize=(8,6))
ax1, ax2 = axs
with np.errstate(over='ignore', invalid='ignore'):
    ax1, _, _ = bifurcation_diagram(ax1, lambda r,x: r-x**2, (1.38,1.41,5000), (-2,2,20), 1000, 10, marker=',', color='k', linestyle='none', alpha=0.1)
ax1.set_xlabel("$r$")
ax1.set_ylabel("$x_n$")
ax1.set_title(r"分叉图和李雅普诺夫指数（红：分叉点，绿：$r_{\infty}$，$r\in[1.38,1.41]$）")
ax1.set_ylim((-0.8,1.5))
ax1.grid()

ax2.plot(r_arr, lambda_r, label=r"$\lambda$")
ax2.grid()
ax2.set_ylim((-2,1))
ax2.set_xlabel("$r$")
ax2.set_ylabel(r"$\lambda$")
ax2.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax1.set_xlim((1.38,1.41))

for ax in axs:
    for rn in res_rn[3:]:
        ax.axvline(rn, linestyle=":", color="r")
    ax.axvline(res_r_inf, linestyle=":", color="g")

display_inline(fig)

可以发现，我们的分叉点和 $r_{\infty}$ 还是比较准确的.

费根鲍姆常数

通过分叉点计算第一费根鲍姆常数

前面我们计算了前几个分叉点位置 $r_n$，而第一费根鲍姆常数定义为：

$$ \delta = \lim_{n\to \infty} \frac{r_{n} - r_{n-1}}{r_{n+1}- r_{n}}$$

我们用前几个分叉点位置来计算，得到的 $\delta$ 随 $n$ 的变化关系为：

delta_exact = 4.669201609
drn = res_rn[1:] - res_rn[:-1]
Y = drn[:-1]
X = drn[1:]
ans = Y/X
n = np.arange(2,len(ans)+2)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(n, ans, "o--r", label=r"$\delta$")
ax.grid()
ax.set_xlabel("n")
ax.set_ylabel(r"$\delta$")
ax.axhline(delta_exact, linestyle=':', color='k',label=f"$\\delta_{{lm}}$={delta_exact:.4f}")
ax.legend()
ax.set_title(r"$\delta$ 收敛图及和 Logisic Map 的 $\delta$ 对比")
display_inline(fig)
display_math(r"\delta\approx",f"{ans[-1]}")

$\displaystyle \delta\approx4.661697250004176$

图中 $\delta_{lm}$ 是课件中 Logistic Map 的费根鲍姆常数，其为：

$$ \delta_{lm} = 4.6692016\dots $$

而我们通过计算得到的映射 $x_{n+1} = r-x^2_n$ 的费根鲍姆常数估计值为：

$$ \delta \approx 4.66169725 $$

两者的相对误差为：

$$ \frac{| \delta-\delta_{lm} |}{\delta_{lm}} \approx 0.16 \% $$

可见二者十分接近！且由于分叉点计算的数量较小，计算出的常数根据收敛趋势判断应该是偏小的，实际的值应该更接近 Logistic Map 的费根鲍姆常数.

通过超稳定环计算第第一、二费根鲍姆常数

前面我们计算费根鲍姆常数是使用了前面找到的分叉点，但是有一个问题限制了精度：临界慢化. 具体而言，在接近分叉点的位置时，$x_n$ 收敛很慢，导致我们需要计算大量步骤后才能得到具体的周期数，同时，也导致了周期数的难以计算和准确性较低，从而限制了找到的分叉点数量，并使得计算的费根鲍姆常数精度一般.

而费根鲍姆常数除了使用分叉点来计算外，还可以使用超稳定环的位置来计算. 具体而言，对于映射 $x_{n+1} = r - x^2_{n}$，其超稳定环出现在 $x^*=0$ 处（此时导数为 $0$），此时费根鲍姆常数定义为：

$$ \delta = \lim_{n\to\infty} \frac{R_n - R_{n-1}}{R_{n+1} -R_n} $$

其中 $R_n$ 为 $2^n$ 的超稳定环出现时的 $r$. 同时还有第二费根鲍姆常数，其定义为：

$$ \alpha = - \lim_{n\to\infty} \frac{d_n}{d_{n+1}} $$

其中 $d_n$ 为超稳定环中离 $x^*=0$ 最近的点到 $x^*=0$ 的 距离.

超稳定环相比分叉点很容易搜索，这是因为我们可以直接以 $0$ 为起点，迭代 $n$ 次，检测最后是否为 $0$（或者在 $0$ 附近），每给定一个 $r$，只需要迭代 $n$ 次，不需要先让 $x_n$ 收敛，从而没有临界慢化问题，大大降低了计算量. 从而使得更精细的搜索成为可能.

而搜索超稳定环还是使用之前的初筛+细化的思路，原先的问题转为寻找 $f^n(0)=0$ 时的 $r$，因此，初筛的思路如下：

• 给定范围 $[a,b]$，其中 $b=a+\Delta r$，$\Delta r$ 为一次搜索的范围
• 计算 $f^n(a)$ 和 $f^n(b)$，如果同号，则移动到下一个步长，如果异号，则记录 $[a,b]$，并进行细化得到 $R_n$
• 新的 $a$ 定义为 $R_n + \Delta r'$，保证已经搜索过的不会再次被搜索.

而细化还是使用二分法，根据 $f^n((a+b)/2)$ 决定区间收缩方向，这里不再赘述.

通过上述方法，我们可以找到各个超稳定环出现的位置：

def refine(func, r_range, tol=1e-15, max_iter=100):
    a, b = r_range[0], r_range[1]
    fa, fb = func(a), func(b)
    for _ in range(max_iter):
        if (b - a) < tol:
            return (a+b)/2
        mid = (a+b)/2
        fmid = func(mid)
        if fa * fmid < 0:
            b = mid
        else:
            a = mid
    else:
        print("到达最大次数后仍未找到.")
    return (a+b)/2

def find_next_r(period, start, end=1.402, step=1e-6):
    def fn(r, n):
        x = 0
        for _ in range(n):
            x = r - x**2
        return x
    eq = lambda r: fn(r, period)

    r = start
    y = eq(r)
    
    guess_range = None
    while r < end:
        r_next = r + step
        if r_next > end: r_next = end
        y_next = eq(r_next)
        if y * y_next < 0:
            guess_range = [r, r_next]
            break
        r = r_next
        y = y_next
        if r == end:
            return None
    return refine(eq, guess_range)

res_Rn = []
res_Rn.extend([0.0])

max_n = 10
for n in range(1, max_n + 1):
    period = 2**n
    start_r = res_Rn[-1] + 1e-9
    r_n = find_next_r(period, start_r)
    if r_n is not None:
        res_Rn.append(r_n)
        print(f"周期 {period:<4} -> R{n+1:<2} = {r_n:.15f}")

周期 2    -> R2  = 1.000000000000000
周期 4    -> R3  = 1.310702641336833
周期 8    -> R4  = 1.381547484432061
周期 16   -> R5  = 1.396945359704560
周期 32   -> R6  = 1.400253081214783
周期 64   -> R7  = 1.400961962944842
周期 128  -> R8  = 1.401113804939775
周期 256  -> R9  = 1.401146325826947
周期 512  -> R10 = 1.401153290849924
周期 1024 -> R11 = 1.401154782546614

对于上面的每一个 $R_n$，我们都可以找到其最后稳定轨道的每一个点，并得到最靠近 $0$ 的点（除了 $0$ 外），从而计算 $d_n$：

@nb.njit()
def calc_dn(r,n,n_warm):
    x = 0
    for _ in range(n_warm):
        x = r - x ** 2
    tmp = np.zeros(n)
    for i in range(n):
        x = r - x**2
        tmp[i] = x
    tmp = np.sort(np.abs(tmp))
    return tmp[1]

res_dn = []

for n, Rn in enumerate(res_Rn):
    period = 2**n
    dn = calc_dn(Rn, period, 1000)
    res_dn.append(dn)
    print(f"周期 {period:<4} -> R{n+1:<2} = {Rn:.15f}, d{n+1:<2} = {dn:.15f}")

周期 1    -> R1  = 0.000000000000000, d1  = 0.000000000000000
周期 2    -> R2  = 1.000000000000000, d2  = 1.000000000000000
周期 4    -> R3  = 1.310702641336833, d3  = 0.407238772670518
周期 8    -> R4  = 1.381547484432061, d4  = 0.163425362208270
周期 16   -> R5  = 1.396945359704560, d5  = 0.065363373186764
周期 32   -> R6  = 1.400253081214783, d6  = 0.026121261182502
周期 64   -> R7  = 1.400961962944842, d7  = 0.010436922716251
周期 128  -> R8  = 1.401113804939775, d8  = 0.004169968060263
周期 256  -> R9  = 1.401146325826947, d9  = 0.001666053648770
周期 512  -> R10 = 1.401153290849924, d10 = 0.000665647584942
周期 1024 -> R11 = 1.401154782546614, d11 = 0.000265949726287

综上，我们可以绘制通过超稳定环得到的两个费根鲍姆常数的收敛过程：

delta_exact = 4.669201609
alpha_exact = -2.502907875
res_Rn = np.array(res_Rn)
dRn = res_Rn[1:] - res_Rn[:-1]
Y1 = dRn[:-1]
X1 = dRn[1:]
ans1 = Y1/X1

res_dn = np.array(res_dn)
ddn = res_dn[1:] - res_dn[:-1]
Y2 = ddn[:-1]
X2 = ddn[1:]
ans2 = -Y2/X2

n = np.arange(2,len(ans1)+2)
fig, ax1 = plt.subplots(figsize=(8,6))
ax2 = ax1.twinx()

ax1.plot(n, ans1, "o--r", label=r"$\delta$")
ax1.grid()
ax1.set_xlabel("n")
ax1.set_ylabel(r"$\delta$")
ax1.axhline(delta_exact, linestyle=':', color='k',label=f"$\\delta_{{lm}}$={delta_exact:.4f}")

ax2.plot(n, ans2, "o--b", label=r"$\delta$")
ax2.set_ylabel(r"$\alpha$")
ax2.axhline(alpha_exact, linestyle=':', color='g',label=f"$\\alpha_{{lm}}$={alpha_exact:.4f}")

lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc='right')

ax1.set_title(r"$\delta,\alpha$ 收敛图及和 Logisic Map 的 $\delta,\alpha$ 对比")
display_inline(fig)
display_math(r"\delta\approx",f"{ans1[-1]}")
display_math(r"\alpha\approx",f"{ans2[-1]}")

$\displaystyle \delta\approx4.669195170819541$

$\displaystyle \alpha\approx-2.5029057378387867$

图中 $\delta_{lm},\alpha_{lm}$ 是课件中 Logistic Map 的两个费根鲍姆常数，其为：

$$ \delta_{lm} = 4.669201609\dots,\quad \alpha_{lm} = -2.502907875$$

而我们通过计算得到的映射 $x_{n+1} = r-x^2_n$ 的费根鲍姆常数估计值为：

$$ \delta \approx 4.66919517, \quad \alpha \approx -2.502905738 $$

两者的相对误差为：

$$ \frac{| \delta-\delta_{lm} |}{\delta_{lm}} \approx 1.4 \times 10^{-6} ,\quad \frac{| \alpha-\alpha_{lm} |}{\alpha_{lm}} \approx 8.5 \times 10^{-7} $$

可见二者及其接近，这进一步说明了映射 $x_{n+1} = r - x^2_{n}$ 的两个费根鲍姆常数和原来的 Logistic Map 相同，同时按照不同方式得到的费根鲍姆常数也相同（通过倍周期分叉点和超稳定环）

选做内容：数值计算李雅普诺夫指数

思路

要数值计算李雅普诺夫指数，就要回归李雅普诺夫指数的原始定义，即定初始点 $x_0$ 和微扰后的点 $x_0 + \delta_0$，分别进行迭代，其距离满足：

$$ \delta_n = \mathrm{e}^{\lambda n} \delta_0 $$

其中 $\lambda$ 即为李雅普诺夫指数.

上式又可以写为：

$$ \lambda = \frac{1}{n} \ln \frac{\delta_n}{\delta_0} $$

而上式成立的条件是 $\delta_n$ 满足（至少近似满足）：

$$ \delta_n = \mathrm{e}^{\lambda n} \delta_0 $$

因此，我们可以先绘制不同情况下 $\delta_n / \delta_0$ 随 $n$ 的变化（取 $r=0.2,r=0.75,r=1.7$），分别对应 $\lambda <0,\lambda=0,\lambda>0$ 三种情况

delta0 = 1e-8
r_arr = np.array([0.2,0.75,1.7])
x0 = np.array([0.1])
x0_p = x0 + delta0
n = 100
dn = np.zeros((n,3))
n_arr = np.arange(n)
for i in range(n):
    x0 = r_arr - x0**2
    x0_p = r_arr - x0_p**2
    dn[i] = np.abs(x0_p - x0)
plt.rcParams['font.family'] = ['DejaVu Sans','SimSun', 'Times New Roman'] # 临时设置字体，消除上标的负号问题
fig, ax = plt.subplots(figsize=(8,6))
for i in range(3):
    ax.plot(n_arr, dn[:,i], ".--", label=f"$r={r_arr[i]:.2f}$")
ax.set_yscale("log")
ax.legend()
ax.grid()
ax.set_xlabel("$n$")
ax.set_ylabel(r"$\delta_n$")
ax.set_title("$\\delta_n$ 随 $n$ 变化，$(\\delta_0 = 10^{-8})$，对数坐标",  fontfamily="SimSun")
plt.rcParams['font.family'] = ['SimSun', 'Times New Roman']
display_inline(fig)

可以看到，对于 $\lambda < 0$ 和 $\lambda=0$，其 $\delta_n$ 确实基本符合 $\delta_n = \mathrm{e}^{\lambda n} \delta_0$（在图中表现为一条直线）

但是，对于 $\lambda >0$，其只在前面的一段距离近似成直线，在一定的 $n$ 后出现了拐点，使得不满足于我们原始的假设（实际上李雅普诺夫指数的定义本身也只在 $\delta$ 是围绕也就是 $\delta \ll 0$ 的情况下成立，当 $\delta_n$ 较大时这个条件就不成立了）

因此，我们的计算思路如下：

• 给定 $x_0$，迭代 $n_{\text{warm}}$ 次，去掉暂态过程.
• 构造 $x_p$ 和 $x'_p = x_p + \delta_0$，迭代计算 $n_{\text{renorm}}$ 次
• 计算此时的 $\delta_n = |x_p - x'_p|$，计算 $\ln(\delta_n / \delta_0)/n_{\text{renorm}}$ 并累加到 $S$ 中，重复第二第三步 $N$ 次
• $\lambda = S/N$

实现和对比

根据上述的思路，我们的实现如下：

@nb.njit()
def calc_lyapunov_1d_by_num(r_arr, n_warm, n_step, n_renorm, x0=0.0, delta0=1e-8):
    f = lambda r,x : r - x**2
    x = x0 * np.ones_like(r_arr)
    for _ in range(n_warm):
        x = f(r_arr, x)
    ncycle = n_step // n_renorm
    sum_lambda = np.zeros_like(r_arr)
    for cycle in range(ncycle):
        x2 = x + delta0
        for _ in range(n_renorm):
            x = f(r_arr, x)
            x2 = f(r_arr, x2)
        delta_n = np.abs(x2-x)
        sum_lambda += np.log(delta_n / delta0) / n_renorm
    return sum_lambda / ncycle

运用上述方法，我们计算全参数范围的李雅普诺夫指数：

r_arr = np.linspace(0,2,10000)
lambda_r = calc_lyapunov_1d_by_num(r_arr, 1000, 200, 10)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, lambda_r, label=r"$\lambda$")
ax.set_ylim((-2,1))
ax.set_xlabel("$r$")
ax.set_ylabel(r"$\lambda$")
ax.set_title("全参数范围的李雅普诺夫指数（数值求解，每 10 次重新缩放）")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

可以看到和原先通过解析方式得到的全参数范围的李雅普诺夫指数相似（除了在 $r=2.0$ 处有一个凸起在原先的图像没有）. 我们将两者叠加显示.

r_arr = np.linspace(0,2,10000)

lambda_r1 = calc_lyapunov_1d_by_num(r_arr, 1000, 200, 10)
lambda_r2 = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            1000,
                            200)
fig, ax = plt.subplots(figsize=(8,6))
ax.set_ylim((-2,1))
ax.set_ylabel(r"$\lambda$")
ax.set_xlabel(r"$r$")
ax.plot(r_arr, lambda_r2, "b", label=r"解析")
ax.plot(r_arr, lambda_r1, "--r", label=r"数值")
ax.set_title("全参数范围的李雅普诺夫指数（不同计算方式对比）")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

我们绘制两者的差的曲线：

fig, ax = plt.subplots(figsize=(8,6))
ax.set_ylabel(r"$\Delta \lambda$")
ax.set_xlabel(r"$r$")
ax.plot(r_arr, lambda_r2-lambda_r1,)
ax.set_title("两种计算方式的差")
ax.set_ylim((-0.5,0.5))
ax.grid()
display_inline(fig)

可以发现，相差较大的部分主要在超稳定环（此时 $\lambda$ 应该为 $-\infty$，但是由于数值精度限制，不可能得到 $-\infty$）.

同时，相差较大的还有在混沌区.（可能是因为在混沌区本身李雅普诺夫指数变化很快？）

而在其他参数，两者均比较吻合.

结论

解析求解周期一和周期二轨道

当 $0 < r < \frac{3}{4}$ 时，系统有稳定的周期一轨道：

$$ \dots \to \frac{\sqrt{4r+1}-1}{2} \to \dots $$

其中吸引域为（Q：可否进一步证明确实是吸引域？需要更好的放缩或者估计？）

$$| x_0 | < \frac{\sqrt{4r+1}+1}{2}$$

当 $\frac{3}{4} < r < \frac{5}{4}$ 时，系统有稳定的周期二轨道：

$$ \dots \to \frac{1+\sqrt{4r-3}}{2} \to \frac{1-\sqrt{4r-3}}{2} \to \dots $$

其中吸引域为

$$| x_0 | < \frac{\sqrt{4r+1}+1}{2}$$

同时，对于周期二，其两条轨道各自的吸引域呈现复杂的结构，其越靠近整体的吸引域边界，吸引域变化越快（Q：可否给出解析通式？）

全周期分叉图

映射 $x_{n+1} = r - x_{n}^2$ 的全周期分叉图如下：

fig, ax = plt.subplots(figsize=(8,6))
with np.errstate(over='ignore', invalid='ignore'):
    ax, r_save, x_save = bifurcation_diagram(ax, lambda r,x: r-x**2, (0,2,5000), (-2,2,20), 1000, 10, marker=',', color='k', linestyle='none', alpha=0.1)
ax.set_xlabel("$r$")
ax.set_ylabel("$x_n$")
ax.set_title(r"分叉图（$x_{n+1} = r-x^2$）")
ax.grid()
display_inline(fig)

其中周期三出现在 $r=\frac{7}{4}$ 处（解析求解），周期五大约出现在 $r \approx 1.625$ 处.

全参数范围李雅普诺夫指数

映射 $x_{n+1} = r - x_{n}^2$ 的全参数范围李雅普诺夫指数如下（通过解析方式求解）：

r_arr = np.linspace(0,2,10000)
lambda_r = calc_lyapunov_1d(lambda r,x: r-x**2,
                            lambda r,x: -2*x,
                            0.1,
                            r_arr,
                            1000,
                            200)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, lambda_r, label=r"$\lambda$")
ax.set_ylim((-2,1))
ax.set_xlabel("$r$")
ax.set_ylabel(r"$\lambda$")
ax.set_title("全参数范围的李雅普诺夫指数")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

倍周期分叉点，$r_{\infty}$，费根鲍姆常数

映射 $x_{n+1} = r - x_n^2$ 的前几个分叉点如下（初筛+0.618 法）：

分叉点	1→2	2→4	4→8	8→16	16→32	32→64
估计值	0.749965	1.249983	1.368091	1.394043	1.399630	1.400828

$r_{\infty}$ 为（初筛+二分法）：

$$ r_{\infty} \approx 1.40115522\dots $$

通过分叉点估计的费根鲍姆常数为：

$$ \delta \approx 4.66169725 $$

与上次作业的 Logistic Map 的费根鲍姆常数的相对偏差为 $0.16\%$

同时，还尝试了通过超稳定环来计算两个费根鲍姆常数，其中前几个超稳定环位置为（初筛+二分法）：

周期	2	4	8	16	32	64	128	256	512	1024
估计值	1.000000	1.310703	1.381547	1.396945	1.400253	1.400962	1.401114	1.401146	1.401153	1.401155

通过超稳定环估计的费根鲍姆常数为：

$$ \delta \approx 4.66919517, \quad \alpha \approx -2.502905738 $$

与 Logistic Map 的相对偏差为：$1.4\times 10^{-6}$ 和 $8.5 \times 10^{-7}$

因此，可以认为映射 $x_{n+1} = r - x_{n}^2$ 的两个费根鲍姆常数和 Logistic Map 相同.

同时，由于分叉点存在临界慢化效应，因此通过分叉点计算费根鲍姆常数相比通过超稳定环计算更麻烦，精度更低，计算量更大.

选作内容：数值计算李雅普诺夫指数

通过数值计算全参数范围的李雅普诺夫指数如下（间隔一定次数重新设 $\delta$）

r_arr = np.linspace(0,2,10000)
lambda_r = calc_lyapunov_1d_by_num(r_arr, 1000, 200, 10)
fig, ax = plt.subplots(figsize=(8,6))
ax.plot(r_arr, lambda_r, label=r"$\lambda$")
ax.set_ylim((-2,1))
ax.set_xlabel("$r$")
ax.set_ylabel(r"$\lambda$")
ax.set_title("全参数范围的李雅普诺夫指数（数值求解，每 10 次重新缩放）")
ax.grid()
ax.axhline(0, linestyle=':', color='k', label=f"$\\lambda=0.0$")
ax.legend()
display_inline(fig)

上面数值计算的结果和解析的相似，除了在 $r=2$ 附近有一个凸起（Why？暂时没想到可能的原因. 数值误差？）

参考资料

• Walter Rudin, 数学分析原理 (Principles of Mathematical Analysis), 3rd ed. §9.3 凝缩原理（也就是压缩映射原理）
• Steven. H. Strogatz, 非线性动力学与混沌 (Nonlinear Dynamics and Chaos), 2nd ed. §10.5, §10:练习 10.6.1
• 部分代码编写时使用 Fitten Code AI 进行补全
• 代码实现中使用 Google Gemini AI 查询了一些函数的用法.

致谢

感谢 9-26 所有在课后讨论李雅普诺夫指数的同学，上述讨论中数值计算李雅普诺夫指数的方法在本次作业中有所体现.（虽然后面的讨论我在符号计算中使用错了函数，将 fn 当成 f 导数实际计算的是 fn 的导数（也就是经过 $n$ 次迭代的函数的导数）导致给出了错误的结论. 在此表示抱歉. QAQ）

附一：问题

（下面的问题除了包含个人能力有限不知道如何解决的问题外，也包含一些受限于时间暂时没有探索的问题）

• 动力学上：
  • 如何分析靠近周期二吸引域边界时出现的不同轨道吸引域的快速变化？这种变化速率是否也满足某种分布（比如给定一个 $r$，其每个“条带”的宽度随条带编号 $k$ 是否也是某种指数关系？）
  • 上面这种轨道吸引域的快速变化是否和后面的混沌有关（相当于混沌区是周期 $n$（$n\to\infty$）的轨道，有 $n$ 条细分的轨道，因此条带宽度非常窄，趋向于 $0$，从而使得由于扰动由于误差等原因很容易在不同条带之间“跳跃”，导致最后行为的随机）
  • $x_{n+1} = r - x_{n}^2$ 和 Logistic Map 有相似分叉图，费根鲍姆常数的原因？一种可能的想法是这两者相当于对 $x$ 进行线性变换 $ax+b$，这个过程是不改变动力学性质的（但是很显然这里 $a,b$ 涉及 $r$）
• 数学上：
  • 对于周期一、周期二整体的吸引域，能否通过合理的估计和放缩，借助压缩映射原理来证明其吸引域？或者有无其他方法证明？是否能证明？
• 数值计算上：
  • 通过数值方式得到的李雅普诺夫指数在 $r=2$ 附近有个凸起的原因？是有具体的动力学原因还是纯粹的数值误差？
  • 如何较好地估计数值计算得到结果的有效位数？（除了看收敛曲线通过不同迭代次数计算结果的差值来估计精度外，是否有其他方式估计？）
  • 如何大致确定每一步计算需要的精度要求？（比如我希望得到的费根鲍姆常数有 $n$ 位有效位数，能否估计要求分叉点有多少位有效位数？再然后如何控制参数使得确实能得到这么多有效位数？）

附二：学习心得

• 解析，数值计算一维系统李雅普诺夫指数的方法，计算一维系统分叉点，超稳定环的方法
• 二分法，0.618 法等求根，优化方法在具体科学计算中的应用
• 不同方法常常能带来不同的计算难度